R version 4.0.3 (2020-10-10) – “Bunny-Wunnies Freak Out”

Packages used for NMDS: vegan (version 2.5-7)

Methods

The document serves as an example of analyses that will be conducted to identify natural variations in benthic communities across Virginia. These NMDS will support the Genus level IBI development process. This analysis is the first run of all of reference sites in Virginia. No West Virginia DEP data is used in this analysis. Several reference sites have been sent back to Virginia DEQ biologists for one more review.

The dataset used is a subset of reference stations collected in Virginia. If stations appeared in the dataset more than 4 times, then the most recent 4 samples were used and the rest removed. Taxa that occurred in the dataset < 10% of the time were removed. The data was log10 +1 transformed. Environmental factors were compiled for each station and used to plot over the NMDS to show environmental variation associated with the community matrix. The envfit function in Vegan was used to plot the continuous environmental variables.

The first step was to read in the reference site bug taxa list and environmental factors dataset for each station. Join the environmental dataset with the bug dataset to account for multiple observations of each station and collection date and time.

Check to make sure the bug and environmental join was successful:

Number of rows in Community Matrix: 854

Number or rows in Environmental Matrix: 858

The data was log10+1 transformed. Rare taxa (<=10%) were removed.

Run NMDS for reference communities

## Run 0 stress 0.1684074 
## Run 1 stress 0.1685779 
## ... Procrustes: rmse 0.003267809  max resid 0.07093681 
## Run 2 stress 0.1684073 
## ... New best solution
## ... Procrustes: rmse 0.0001266153  max resid 0.002595457 
## ... Similar to previous best
## Run 3 stress 0.1694493 
## Run 4 stress 0.1687056 
## ... Procrustes: rmse 0.00336004  max resid 0.07511734 
## Run 5 stress 0.1685908 
## ... Procrustes: rmse 0.002140328  max resid 0.06192395 
## Run 6 stress 0.1686516 
## ... Procrustes: rmse 0.002420567  max resid 0.06191808 
## Run 7 stress 0.1684666 
## ... Procrustes: rmse 0.001172967  max resid 0.0330259 
## Run 8 stress 0.1685785 
## ... Procrustes: rmse 0.003736797  max resid 0.07435475 
## Run 9 stress 0.1687616 
## ... Procrustes: rmse 0.004266263  max resid 0.07264733 
## Run 10 stress 0.1685902 
## ... Procrustes: rmse 0.002139672  max resid 0.06194481 
## Run 11 stress 0.1687013 
## ... Procrustes: rmse 0.004134911  max resid 0.07368102 
## Run 12 stress 0.1683998 
## ... New best solution
## ... Procrustes: rmse 0.00245399  max resid 0.07099414 
## Run 13 stress 0.1686374 
## ... Procrustes: rmse 0.002434842  max resid 0.06176666 
## Run 14 stress 0.168519 
## ... Procrustes: rmse 0.002543314  max resid 0.0737268 
## Run 15 stress 0.1685776 
## ... Procrustes: rmse 0.002132679  max resid 0.06176917 
## Run 16 stress 0.1685912 
## ... Procrustes: rmse 0.003247928  max resid 0.07120008 
## Run 17 stress 0.1685785 
## ... Procrustes: rmse 0.003628517  max resid 0.07125442 
## Run 18 stress 0.1684073 
## ... Procrustes: rmse 0.00245374  max resid 0.07099415 
## Run 19 stress 0.1685178 
## ... Procrustes: rmse 0.00361  max resid 0.0767033 
## Run 20 stress 0.1685776 
## ... Procrustes: rmse 0.002137212  max resid 0.06179338 
## Run 21 stress 0.1695664 
## Run 22 stress 0.1694941 
## Run 23 stress 0.1687015 
## ... Procrustes: rmse 0.003358745  max resid 0.07518244 
## Run 24 stress 0.1685774 
## ... Procrustes: rmse 0.003783563  max resid 0.07639676 
## Run 25 stress 0.1696084 
## Run 26 stress 0.1685773 
## ... Procrustes: rmse 0.00213494  max resid 0.06176447 
## Run 27 stress 0.1684077 
## ... Procrustes: rmse 0.002456996  max resid 0.0710458 
## Run 28 stress 0.1684075 
## ... Procrustes: rmse 0.002453015  max resid 0.07097477 
## Run 29 stress 0.1696052 
## Run 30 stress 0.168407 
## ... Procrustes: rmse 0.002451901  max resid 0.07099345 
## Run 31 stress 0.168637 
## ... Procrustes: rmse 0.002428719  max resid 0.06177029 
## Run 32 stress 0.1687438 
## ... Procrustes: rmse 0.003274938  max resid 0.07480915 
## Run 33 stress 0.1694864 
## Run 34 stress 0.1684581 
## ... Procrustes: rmse 0.001195828  max resid 0.03324885 
## Run 35 stress 0.1686511 
## ... Procrustes: rmse 0.003441384  max resid 0.07110878 
## Run 36 stress 0.1687412 
## ... Procrustes: rmse 0.004113679  max resid 0.07639219 
## Run 37 stress 0.1684821 
## ... Procrustes: rmse 0.00275688  max resid 0.07083849 
## Run 38 stress 0.1685773 
## ... Procrustes: rmse 0.00376374  max resid 0.07567305 
## Run 39 stress 0.1685907 
## ... Procrustes: rmse 0.003244964  max resid 0.07116412 
## Run 40 stress 0.1687676 
## ... Procrustes: rmse 0.004319747  max resid 0.07585751 
## Run 41 stress 0.1684656 
## ... Procrustes: rmse 0.002715933  max resid 0.07085339 
## Run 42 stress 0.1685902 
## ... Procrustes: rmse 0.003242766  max resid 0.07115944 
## Run 43 stress 0.1686395 
## ... Procrustes: rmse 0.002090333  max resid 0.04484982 
## Run 44 stress 0.1687612 
## ... Procrustes: rmse 0.003495713  max resid 0.07337875 
## Run 45 stress 0.1687013 
## ... Procrustes: rmse 0.003326373  max resid 0.07395997 
## Run 46 stress 0.1685778 
## ... Procrustes: rmse 0.00213211  max resid 0.06177182 
## Run 47 stress 0.1684074 
## ... Procrustes: rmse 0.002451144  max resid 0.07098344 
## Run 48 stress 0.1684001 
## ... Procrustes: rmse 0.0001359149  max resid 0.002198067 
## ... Similar to previous best
## *** Solution reached
## 
## Call:
## metaMDS(comm = NMDSten[, 6:83], k = 3, trymax = 1000) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     NMDSten[, 6:83] 
## Distance: bray 
## 
## Dimensions: 3 
## Stress:     0.1683998 
## Stress type 1, weak ties
## Two convergent solutions found after 48 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'NMDSten[, 6:83]'

Envfit results from Vegan package :

## 
## ***VECTORS
## 
##                     NMDS1    NMDS2     r2 Pr(>r)   
## JulianDate       -0.54082  0.84114 0.0745  0.658   
## Latitude          0.87395 -0.48602 0.2561  0.192   
## Longitude         0.98919  0.14663 0.2643  0.189   
## US_L3CODE        -0.13768 -0.99048 0.2884  0.146   
## Order            -0.96417 -0.26530 0.7460  0.002 **
## EDASOrder        -0.96417 -0.26530 0.7460  0.002 **
## totalArea_sqMile -0.72068 -0.69327 0.4557  0.049 * 
## ELEVMEAN         -0.89598 -0.44410 0.5281  0.013 * 
## SLPMEAN          -0.42780 -0.90387 0.7659  0.002 **
## wshdRain_mmyr    -0.70382 -0.71037 0.4590  0.049 * 
## siteRain_mmyr    -0.65103  0.75905 0.2733  0.182   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                                                          NMDS1   NMDS2
## SeasonFall                                             -0.5647  0.6258
## SeasonOutside Sample Window                            -0.5792  0.1746
## SeasonSpring                                           -0.5968  0.3210
## GradientMACS                                           -0.5818  0.4412
## CoastalTRUE                                            -0.5818  0.4412
## US_L3NAMEMiddle Atlantic Coastal Plain                 -0.5219  0.6919
## US_L3NAMESoutheastern Plains                           -0.6151  0.3019
## US_L4NAMEChesapeake-Pamlico Lowlands and Tidal Marshes -0.5791  0.3916
## US_L4NAMEMid-Atlantic Flatwoods                        -0.4837  0.8921
## US_L4NAMERolling Coastal Plain                         -0.6151  0.3019
## ASSESS_REGPRO                                          -0.6086  0.3182
## ASSESS_REGTRO                                          -0.4837  0.8921
## TidalN                                                 -0.5823  0.4494
## TidalT                                                 -0.5791  0.3916
## VAHUSBCB                                               -0.5022  0.5568
## VAHUSBCU                                               -0.6093  0.5686
## VAHUSBYO                                               -0.6202  0.1344
## BasinChes. Bay and Small Coastal Basin                 -0.5022  0.5568
## BasinChowan and Dismal Swamp River Basin               -0.6093  0.5686
## BasinYork River Basin                                  -0.6202  0.1344
## Basin_CodeChowan-Dismal                                -0.6093  0.5686
## Basin_CodeSmall Coastal                                -0.5022  0.5568
## Basin_CodeYork                                         -0.6202  0.1344
## CountyCityNameCaroline                                 -0.6891  0.0814
## CountyCityNameDinwiddie                                -0.7258  0.4332
## CountyCityNameHanover                                  -0.5513  0.1875
## CountyCityNameKing and Queen                           -0.4252  0.7219
## CountyCityNameNorthumberland                           -0.5791  0.3916
## CountyCityNameSouthampton                              -0.4837  0.8921
## CountyCityNameSussex                                   -0.7531 -0.1309
## WQS_CLASSIII                                           -0.6579  0.1630
## WQS_CLASSVII                                           -0.5395  0.5957
## WQS_SPSTDS                                             -0.5769  0.3975
## WQS_SPSTDSNEW-21                                       -0.6455  1.0095
## WQS_PWS                                                -0.5818  0.4412
## WQS_TROUT                                              -0.5818  0.4412
## WQS_TIER_III                                           -0.5818  0.4412
## 
## Goodness of fit:
##                    r2 Pr(>r)   
## Season         0.1958  0.341   
## Gradient       0.0000  1.000   
## Coastal        0.0000  1.000   
## US_L3NAME      0.2663  0.075 . 
## US_L4NAME      0.4267  0.043 * 
## ASSESS_REG     0.4189  0.017 * 
## Tidal          0.0030  0.978   
## VAHUSB         0.2901  0.148   
## Basin          0.2901  0.148   
## Basin_Code     0.2901  0.148   
## CountyCityName 0.8557  0.010 **
## WQS_CLASS      0.3332  0.034 * 
## WQS_SPSTDS     0.1814  0.133   
## WQS_PWS        0.0000  1.000   
## WQS_TROUT      0.0000  1.000   
## WQS_TIER_III   0.0000  1.000   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## 840 observations deleted due to missingness

Plot with Station IDs

Plot with Axis 3

Plot with Species

Season

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$Season,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       Fall   Outside Sample Window Spring
## delta 0.6419 0.6564                0.6475
## n     414    13                    427   
## 
## Chance corrected within-group agreement A: 0.03778 
## Based on observed delta 0.6449 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

Season with continuous variable vectors

Ecoregion

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$US_L3NAME,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       Blue Ridge Central Appalachians Middle Atlantic Coastal Plain
## delta 0.5768     0.5884               0.507                        
## n     154        40                   25                           
##       Northern Piedmont Piedmont Ridge and Valley Southeastern Plains
## delta 0.6265            0.5919   0.6382           0.5931             
## n     101               138      259              137                
## 
## Chance corrected within-group agreement A: 0.09754 
## Based on observed delta 0.6049 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

Basin

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$Basin_Code,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       Appomattox Chowan-Dismal James-Lower James-Middle James-Upper New   
## delta 0.5512     0.6223        0.5695      0.5932       0.6337      0.6312
## n     15         33            36          61           120         78    
##       Potomac-Lower Potomac-Shenandoah Rappahannock Roanoke Small Coastal
## delta 0.6619        0.6423             0.648        0.5886  0.528        
## n     31            42                 144          87      21           
##       Tennessee-Big Sandy Tennessee-Clinch Tennessee-Holston Yadkin York  
## delta 0.5827              0.6114           0.6064            0.5361 0.6313
## n     16                  50               53                4      63    
## 
## Chance corrected within-group agreement A: 0.07782 
## Based on observed delta 0.6181 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

Admin Region

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$ASSESS_REG,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       BRRO   NRO    PRO    SWRO   TRO    VRO   
## delta 0.6224 0.6706 0.5961 0.6149 0.4894 0.6291
## n     250    195    122    156    21     110   
## 
## Chance corrected within-group agreement A: 0.06618 
## Based on observed delta 0.6259 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

Sample Method

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$Gradient,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       MACS   Riffle
## delta 0.5978 0.637 
## n     163    691   
## 
## Chance corrected within-group agreement A: 0.06075 
## Based on observed delta 0.6295 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

BioRegion- Not in dataset

Stream Order

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$Order,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       1      2      3      4      5     6     
## delta 0.6472 0.6598 0.6353 0.6293 0.592 0.4265
## n     277    248    172    114    41    2     
## 
## Chance corrected within-group agreement A: 0.0408 
## Based on observed delta 0.6429 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

Water Quality Standard Class

## 
## Call:
## mrpp(dat = bugsnms[, 6:83], grouping = samplescoresenv$WQS_CLASS,      distance = "bray") 
## 
## Dissimilarity index: bray 
## Weights for groups:  n 
## 
## Class means and counts:
## 
##       III    IV     V      VI     VII   
## delta 0.6603 0.6257 0.6159 0.5843 0.5843
## n     331    185    87     209    42    
## 
## Chance corrected within-group agreement A: 0.06609 
## Based on observed delta 0.6259 and expected delta 0.6702 
## 
## Significance of delta: 0.001 
## Permutation: free
## Number of permutations: 999

Natural Trout Water WQS?